# Temas para os gráficos
theme.base <- theme_minimal(base_size = 11) +
theme(
axis.text = element_text(size = 8),
plot.title = element_text(hjust = 0.5, size = 11, face = "bold"),
plot.subtitle = element_text(hjust = 0.5, size = 9),
plot.caption = element_text(size = 8),
axis.title = element_text(size = 8),
legend.title = element_text(size = 8),
panel.grid.major = element_line(colour = "grey90", linewidth = 0.5),
panel.grid.minor = element_line(
colour = adjustcolor("grey90", alpha.f = 0.4),
linewidth = 0.5
),
panel.border = element_blank(),
panel.background = element_blank(),
plot.background = element_blank(),
axis.line.x = element_line(colour = "grey"),
axis.line.y = element_line(colour = "grey"),
)
theme.no_legend <- theme(legend.position = "none")
theme.no_grid <- theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
theme.no_axis <- theme(axis.line.x = element_blank(), axis.line.y = element_blank())
theme.no_axis_title <- theme(axis.title.x = element_blank(), axis.title.y = element_blank())
# Theme for timeseries
theme.ts <- theme.base +
theme(
axis.text.x = element_text(angle = 45, hjust = 1),
panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank()
) +
theme.no_legend
plotly.base <- function(p) {
p %>%
layout(margin = list(b = 60, t = 80)) %>%
config(mathjax = 'cdn')
}# Função para salvar e carregar resultados
cache_dados <- function(chave, funcao_geradora) {
c <- paste("cache", chave, sep = "/")
c <- paste(c, "rds", sep = ".")
cache <- file.exists(c)
if (!cache) {
resultado <- funcao_geradora()
saveRDS(resultado, c)
} else {
resultado <- readRDS(c)
}
resultado
}O βAR(1) é um modelo de séries temporais autorregressiva de ordem 1 com distribuição Beta – onde a variável de interesse está restrita no intervalo (0, 1) –, este modelo foi proposto por Rocha e Cribari-Neto:
Sejam a série temporal \(Y_{t} = \left\{y_{1}, y_{2}, \ldots, y_{t} | y_{n} \in (0, 1)\right\}\), o conjunto \(l\)-dimensional de covariáveis \(\mathbf{X_{t}}\), e o vetor de parâmetros \(\beta = \left\{\beta_{1}, \beta_{2}, \ldots, \beta_{l}\right\}\), temos o modelo βAR(1) é definido como:
\[ g(\mu_t) = \alpha + \mathbf{X_t'}\beta + \phi \left(g(Y_{t-1}) - \mathbf{X_{t-1}\beta}\right) \]
Onde \(g: (0,1) \rightarrow \mathbb{R}\) é a função de ligação, \(\alpha\) é o intercepto, \(\phi\) é o parâmetro de autorregressão, e \(\mu_t\) é a média condicional da distribuição Beta. A função de ligação mais usual é a função logit.
Neste trabalho, objetiva-se aplicar métodos de controle de processos βAR(1) para detectar mudanças no processo. Para isso, utilizamos que, quando modelada corretamente, os resíduos de uma série temporal, como a βAR(1), são independentes e normalmente distribuídos com média zero e variância constante.
Desta forma, quando o processo sofre uma mudança, espera-se que os resíduos não sejam mais independentes e que a média e a variância dos resíduos sejam diferentes. Assim, podemos utilizar métodos de controle de processos para detectar essas mudanças.
Olhando sob a perspectiva de teste de hipóteses, temos:
Além disso, utilizaremos um nível de significância de 0.05 para os testes de hipóteses, o que nos dá um quantil de 1.96 para a distribuição normal padrão.
A simulação dos dados será feita a partir da biblioteca BTSR: Bounded Time Series Regression.
nH0 <- 100
nH1 <- 200
phi_parametro <- 0.2
alpha_teste <- 0.05
quantil_teste <- qnorm(1 - alpha_teste / 2)
coeficientes <- function(phi) {
# βARMA: o modelo de Rocha e Cribari-Neto (2009, 2017)
# é obtido definindo `coefs$d = 0`
# e `d = FALSE` e `error.scale = 1` (escala preditiva)
#
# ν (nu): parâmetro de precisão, quanto maior o seu valor,
# menor a variância condicional (para μ_t fixo).
# O pacote BTSR define como padrão ν = 20.
list(
alpha = 0,
nu = 20,
phi = phi,
d = 0,
nu = 20
)
}
barma.sim <- function(n,
phi,
seed = NULL,
y.start = NULL) {
BARFIMA.sim(
n = n,
coefs = coeficientes(phi),
y.start = y.start,
error.scale = 1,
complete = F,
seed = ifelse(is.null(seed), sample(1:1000, 1), seed)
)
}
barma.phi_estimado <- function(yt,
alpha = 0,
nu = 20,
phi = 0.1) {
BARFIMA.fit(
yt = yt,
start = list(alpha = alpha, nu = nu, phi = phi),
p = 1,
# Odem do polinômio AR
d = FALSE,
error.scale = 1,
report = F
)$coefficients["phi"][[1]]
}
barma.residuos <- function(yt, phi_estimado) {
BARFIMA.extract(
yt = yt,
coefs = coeficientes(phi_estimado),
llk = F,
info = F,
error.scale = 1
)$error
}# Função para forçar a execução de `BARFIMA.extract` até que ela funcione
so_quero_que_funcione <- function(func,
debug = F,
max = 5,
...) {
FINALMENTE <- F
contador <- 0
while (!FINALMENTE) {
contador <- contador + 1
if (debug)
print(paste("Tentativa:", contador))
resultado <- tryCatch({
func(...)
}, error = function(err) {
if (contador >= max) {
stop("Deu ruim, número máximo de tentativas atingido")
}
if (any(grepl("\\.btsr\\.extract\\(", deparse(err$trace$call)))) {
# Ignora erro de extração de resíduos
return(NULL)
}
stop(err)
})
if (!is.null(resultado)) {
FINALMENTE <- T
}
}
return(resultado)
}# Função auxiliar para gerar os dados
# Exemplo 1:
# ```r
# teste1 <- gerador_monte_carlo(
# parametros = list(
# lambda = seq(0.1, 0.4, by = 0.1)
# ),
# numero_de_execucoes = 2
# )
#
# teste1
# ```
#
# Exemplo 2:
# ```r
# teste2 <- gerador_monte_carlo(
# parametros = expand.grid(
# desvio_detectavel = seq(0.6, 1.8, by = 0.2),
# intervalo_de_decisao = seq(4, 6, by = 1)
# )
# numero_de_execucoes = 2
# )
#
# teste2
# ```
#
# Total de execuções é:
# `execucoes = numero_de_execucoes * length(novos_phis) * nrow(parametros)`
#
gerador_monte_carlo <- function(parametros = NULL,
numero_de_execucoes = 200) {
# Verifica se `parametros` é um data frame, caso contrário converte para um data frame
if (!is.null(parametros) && !is.data.frame(parametros)) {
parametros <- as.data.frame(parametros)
}
novos_phis <- seq(0.2, 0.6, by = 0.1)
# Expande a grade de parâmetros para cobrir todas as combinações
result <- expand.grid(k = 1:numero_de_execucoes, h1_phi = novos_phis)
if (!is.null(parametros)) {
result <- merge(result, parametros, all = TRUE)
}
result %>%
mutate(id = row_number()) %>%
rowwise() %>%
mutate(
# H0
## Gera amostra de controle
h0_amostras = list(barma.sim(
n = nH0, phi = phi_parametro, seed = id
)),
## Estima o valor de phi da amostra de controle
h0_phi = barma.phi_estimado(h0_amostras),
## Calcula os resíduos da amostra de controle
h0_residuos = list(barma.residuos(h0_amostras, h0_phi)),
# H1
## Gera a amostra subsequente
h1_controle = list(
barma.sim(
n = nH1,
phi = h1_phi,
# última observação da amostra de controle
y.start = h0_amostras[nH0],
seed = id + 1337E4
)
),
## Calcula os resíduos da amostra subsequente
h1_residuos = list(so_quero_que_funcione(\() barma.residuos(
h1_controle, h0_phi
)))
)
}Vamos verificar nossa implementação do gerador de Monte Carlo.
teste_montecarlo <- cache_dados("teste-montecarlo", function() {
gerador_monte_carlo(parametros = list(nH0 = c(
rep(25, 20), rep(50, 20), rep(100, 20), rep(200, 20)
)),
numero_de_execucoes = 1) %>%
select(nH0, h0_phi, h1_phi) %>%
mutate(# Calcula a diferença entre os valores de phi
# `phi_parametro`: valor de phi definido para a amostra de controle
# `h0_phi`: valor de phi estimado para a amostra de controle. Espera-se que seja igual a `phi_parametro`
diferenca = h0_phi - phi_parametro)
})teste_montecarlo %>%
group_by(nH0) %>%
summarise(
simulações = n(),
mean = mean(diferenca),
var = var(diferenca),
min = min(diferenca),
max = max(diferenca),
.groups = "drop"
)ggplotly(
teste_montecarlo %>%
ggplot(aes(x = factor(nH0), y = diferenca)) +
geom_boxplot() +
geom_hline(
yintercept = 0,
color = "red",
linetype = "dashed"
) +
annotate(
geom = "text",
x = 0.55,
y = 0 + 0.01,
label = "0.00",
color = "red"
) +
labs(
x = "Número de execuções",
y = "Diferença entre Φ e Φ₀",
title = paste0("Diferença entre Φ e Φ₀")
) +
theme.base
) %>%
plotly.baseO EWMA (Exponential Weighted Moving Average) é um método de controle de processos que utiliza uma média móvel ponderada exponencialmente para detectar mudanças no processo. Sendo definido por:
\[ z_i = \lambda y_i + (1 - \lambda) z_{i-1} \]
Onde, \(z_i\) é a estatística de controle no instante \(i\), \(y_i\) é a observação no instante \(i\), \(\lambda\) é o fator de suavização e \(z_{i-1}\) é a estatística de controle no instante anterior.
O fator \(\lambda\) é uma constante definida no intervalo \((0, 1]\) e seu valor inicial (para \(i = 1\)) é definido como a média do processo, tal que \(z_0 = \mu_0\).
A estatística de controle, \(z_i\), é comparada com os limites de controle, \(\text{LCS}\) e \(\text{LCI}\), sendo é definido como:
\[ \text{LC} = \mu_0 \pm L \sigma \sqrt{\frac{\lambda}{2 - \lambda} \left[1 - \left(1 - \lambda\right)^{2i}\right]} \]
Aqui, utilizamos o pacote qcc para implementar o
EWMA.
ewma_qcc <- function(amostra_inicial,
lambda,
amostra_subsequente,
nsigmas = 1.96) {
ew <- ewma(
amostra_inicial,
lambda = lambda,
nsigmas = nsigmas,
plot = F,
newdata = amostra_subsequente
)
registros <- (nH0 + 1):(nH0 + nH1)
ewma <- as.numeric(ew$y[registros])
LI <- ew$limits[, 1][registros]
LS <- ew$limits[, 2][registros]
fora_de_controle <- ewma < LI | ewma > LS
total_fora_de_controle <- sum(fora_de_controle, na.rm = TRUE)
fracao_fora_de_controle <- total_fora_de_controle / length(ewma)
list(
ewma = ewma,
LI = LI,
LS = LS,
fora_de_controle = fora_de_controle,
total_fora_de_controle = total_fora_de_controle,
fracao_fora_de_controle = fracao_fora_de_controle
)
}α = 0.05Vamos analisar buscar o melhor valor de \(\lambda\) para o EWMA, mantendo um nível de significância constante, \(\alpha = 0.05\).
ewma_monte_carlo <- cache_dados("simulacao-ewma", function() {
gerador_monte_carlo(parametros = list(lambda = seq(0.1, 0.9, by = 0.1))) %>%
mutate(
qcc = list(ewma_qcc(h0_residuos, lambda, h1_residuos)),
ewma = list(qcc$ewma),
LI = list(qcc$LI),
LS = list(qcc$LS),
fora_de_controle = list(qcc$fora_de_controle),
total_fora_de_controle = qcc$total_fora_de_controle,
fracao_fora_de_controle = qcc$fracao_fora_de_controle
) %>%
select(-qcc)
})
ewma_monte_carlo$lambda <- as.factor(ewma_monte_carlo$lambda)λ = 0.2Vamos analisar o EWMA para \(\lambda = 0.2\). Aqui, vamos considerar apenas a primeira execução.
ggplotly(
ewma_monte_carlo %>%
filter(k == 1 & lambda == 0.2) %>% # Apenas a primeira execução
select(h1_phi, ewma, LI, LS, fora_de_controle, lambda) %>%
rename(`Φ₁` = h1_phi) %>%
mutate(observacao = list(1:nH1)) %>%
unnest(cols = c(
`Φ₁`, observacao, ewma, LI, LS, fora_de_controle
)) %>%
ggplot() +
geom_line(aes(
x = observacao, y = LI, color = "Limite Inferior"
), linetype = "dashed") +
geom_line(aes(
x = observacao, y = LS, color = "Limite Superior"
), linetype = "dashed") +
geom_line(aes(
x = observacao, y = ewma, color = "EWMA"
)) +
geom_point(
data = . %>% filter(fora_de_controle),
aes(x = observacao, y = ewma, color = "Fora de controle"),
size = 1,
shape = 4
) +
labs(x = "Observação", y = "Caracterísitca", color = "Legenda") +
scale_color_manual(
values = c(
"Limite Inferior" = adjustcolor("red", alpha.f = 0.5),
"Limite Superior" = adjustcolor("red", alpha.f = 0.5),
"EWMA" = adjustcolor("black", alpha.f = 0.8),
"Fora de controle" = adjustcolor("#f42f2f", alpha.f = 0.8)
)
) +
labs(title = "EWMA para resíduos βAR(1), com λ = 0.2", ) +
theme.base +
theme(legend.position = "bottom", legend.title = element_blank()) +
facet_wrap(vars(`Φ₁`), ncol = 2, labeller = "label_both")
) %>%
plotly.baseResumo da fração de pontos fora de controle para diferentes valores de \(\lambda\) e \(\Phi_1\).
ewma_monte_carlo_resumo <- ewma_monte_carlo %>%
group_by(lambda, h1_phi) %>%
summarise(
mean = mean(fracao_fora_de_controle),
min = min(fracao_fora_de_controle),
max = max(fracao_fora_de_controle),
.groups = "drop"
)
datatable(
ewma_monte_carlo_resumo %>%
mutate(across(where(is.numeric), \(x) round(x, digits = 4))),
caption = "Pontos fora de controle por Φ₁ e λ",
colnames = c("Valores de λ", "Valores de Φ₁", "Média", "Mínimo", "Máximo")
)Pela análise gráfica das médias de FPFCs (Fração de Pontos Fora de Controle) notamos que \(\lambda = 0.7\) é o melhor valor para detectar mudanças no processo.
ggplotly(
ewma_monte_carlo_resumo %>%
ggplot(aes(
x = h1_phi, y = mean, color = lambda
)) +
geom_line() +
geom_point() +
labs(
x = "Valores de Φ",
y = "Fração de pontos fora de controle",
color = "Valores de λ",
title = "FPFCs por Φ e λ"
) +
theme.base
) %>%
plotly.baseE, a distribuição dos PFCs.
ggplotly(
ewma_monte_carlo %>%
ggplot(aes(
x = factor(h1_phi),
y = fracao_fora_de_controle,
fill = lambda
)) +
geom_boxplot() +
labs(
x = "Valores de Φ",
y = "Fração de pontos fora de controle",
fill = "Valores de λ",
title = "FPFCs por Φ e λ"
) +
theme.base
) %>%
plotly.base %>%
layout(boxmode = 'group')Segundo Montgomery (2009), o EWMA-AR é uma extensão do EWMA que utiliza um modelo autorregressivo para prever a próxima observação.
Assim, temos que \(\lambda \in (0, 1]\), sendo que, a previsão para a observação \(x_{t+1}\) é dada por \(\hat{x}_{t+1}(t)=z_{t} = \lambda x_t + (1 - \lambda) z_{t-1}\).
De acordo com Montgomery (2009), podemos encontrar um valor ótimo para \(\lambda\) através da minimização da soma dos quadrados dos resíduos.
E, ainda, temos que os erros de previsão são dados por \(e_t = x_t - \hat{x}_t(t-1)\) conforme a Eq. 10.16 (Montgomery, 2009).
Ou, de outra forma
\[ \hat{\lambda} = \min \left( \underset{\lambda\in(0,1]}{\arg\min}\,\text{Err}(\lambda) \right) \]
onde \(\hat{\lambda}\) é o \(\lambda\) ótimo por simulação, e \(\text{Err}(\lambda) = \sum_{t=1}^{n} e_t^2\).
Segundo, também Montgomery (2009), podemos estimar o valor de \(\sigma^2\) para o modelo EWMA-AR como \(\sigma^2 = \frac{\sum (\text{err}_i^2|_\lambda)}{n}\). Onde \(\text{err}|_\lambda\) são os resíduos do modelo EWMA-AR para o melhor valor de \(\lambda\) encontrado.
# Função para computar resíduos do EWMAAR
ewma_ar_residuos <- function(dados, ewma) {
n <- length(dados)
residuos <- numeric(n)
residuos[1] <- dados[1]
for (i in 2:n) {
residuos[i] <- dados[i] - ewma[i - 1]
}
return(residuos)
}
# Função para computar o EWMA-AR
ewma_ar <- function(dados,
lambda,
x0 = NULL,
desvio = NULL) {
desv <- ifelse(is.null(desvio), sqrt(sd(dados)), desvio)
n <- length(dados)
final <- ifelse(is.null(x0), n, n + 1)
ewma_serie <- numeric(final)
ewma_serie[1] <- ifelse(is.null(x0), dados[1], x0)
for (i in 2:final) {
ewma_serie[i] <- lambda * dados[i] + (1 - lambda) * ewma_serie[i - 1]
}
ewma <- ewma_serie[ifelse(is.null(x0), 1, 2):final]
LI <- ewma - 1.65 * desv
LS <- ewma + 1.65 * desv
fora_de_controle <- dados < LI | dados > LS
total_fora_de_controle <- sum(fora_de_controle, na.rm = TRUE)
fracao_fora_de_controle <- total_fora_de_controle / length(ewma)
return(
list(
ewma = ewma,
residuos = ewma_ar_residuos(dados, ewma_serie),
LI = LI,
LS = LS,
fora_de_controle = fora_de_controle,
total_fora_de_controle = total_fora_de_controle,
fracao_fora_de_controle = fracao_fora_de_controle
)
)
}lambda_otimo_optimize <- function(amostra_inicial) {
lambda <- NA
minimo <- Inf
for (x in c(seq(0.01, 0.1, by = 0.01), seq(0.1, 0.9, by = 0.1))) {
ew <- ewma_ar(amostra_inicial, lambda = x)
erros <- ew$residuos
soma_quadrados <- sum(erros ^ 2)
if (soma_quadrados < minimo) {
minimo <- soma_quadrados
lambda <- x
}
}
return(list(lambda = lambda, minimo = minimo))
}
ewmar_lambda_otimo <- cache_dados("ewmaar-lambda-otimo", \() {
data.frame(id = 1:100) %>%
rowwise() %>%
mutate(
h0_amostra = list(barma.sim(n = nH0, phi = phi_parametro)),
h0_phi = list(barma.phi_estimado(h0_amostra)),
h0_residuos = list(barma.residuos(h0_amostra, h0_phi)),
ewma_ar = list(lambda_otimo_optimize(h0_residuos)),
lambda = ewma_ar$lambda,
minimo = ewma_ar$minimo
) %>%
select(-ewma_ar)
})ewmar_lambda_otimo %>%
group_by() %>%
summarise(
mean_lambda = mean(lambda),
mean_err = mean(minimo),
var_err = var(minimo),
min_err = min(minimo),
max_err = max(minimo),
.groups = "drop"
)Por simulação encontramos um \(\lambda\) ótimo de \(\sim 0.07\). Com erro médio de \(\sim 22\), o que nos dá um desvio de \(22 / 100 \simeq 0.2\).
ewmaar_monte_carlo <- cache_dados("simulacao-ewma-ar", function() {
gerador_monte_carlo(parametros = list(lambda = seq(0.01, 0.1, by = 0.01))) %>%
mutate(
ewma_ar = list(ewma_ar(
h1_controle, lambda, h0_amostras[nH0], sd(h0_amostras)
)),
ewma = list(ewma_ar$ewma),
LI = list(ewma_ar$LI),
LS = list(ewma_ar$LS),
fora_de_controle = list(ewma_ar$fora_de_controle),
total_fora_de_controle = ewma_ar$total_fora_de_controle,
fracao_fora_de_controle = ewma_ar$fracao_fora_de_controle
) %>%
select(-ewma_ar)
})
ewmaar_monte_carlo$lambda <- as.factor(ewmaar_monte_carlo$lambda)λ = 0.07Vamos analisar o EWMA-AR para \(\lambda = 0.07\). O resultado encontrado anteriormente.
ggplotly(
ewmaar_monte_carlo %>%
filter(k == 1 & lambda == 0.07) %>% # Apenas a primeira execução
select(h1_phi, h1_controle, LI, LS, fora_de_controle, lambda) %>%
rename(`Φ₁` = h1_phi) %>%
mutate(observacao = list(1:nH1)) %>%
unnest(cols = c(
`Φ₁`, observacao, h1_controle, LI, LS, fora_de_controle
)) %>%
ggplot() +
geom_line(aes(
x = observacao, y = LI, color = "Limite Inferior"
), linetype = "dashed") +
geom_line(aes(
x = observacao, y = LS, color = "Limite Superior"
), linetype = "dashed") +
geom_line(aes(
x = observacao, y = h1_controle, color = "Série"
)) +
geom_point(
data = . %>% filter(fora_de_controle),
aes(x = observacao, y = h1_controle, color = "Fora de controle"),
size = 1,
shape = 4
) +
labs(x = "Observação", y = "Caracterísitca", color = "Legenda") +
scale_color_manual(
values = c(
"Limite Inferior" = adjustcolor("red", alpha.f = 0.5),
"Limite Superior" = adjustcolor("red", alpha.f = 0.5),
"Série" = adjustcolor("black", alpha.f = 0.8),
"Fora de controle" = adjustcolor("#f42f2f", alpha.f = 0.6)
)
) +
labs(title = "EWMA para resíduos βAR(1), com λ = 0.08", ) +
theme.base +
theme(legend.position = "bottom", legend.title = element_blank()) +
facet_wrap(vars(`Φ₁`), ncol = 2, labeller = "label_both")
) %>%
plotly.baseResumo da fração de pontos fora de controle para diferentes valores de \(\lambda\) e \(\Phi_1\).
ewmaar_monte_carlo_resumo <- ewmaar_monte_carlo %>%
group_by(lambda, h1_phi) %>%
summarise(
mean = mean(fracao_fora_de_controle),
min = min(fracao_fora_de_controle),
max = max(fracao_fora_de_controle),
.groups = "drop"
)
datatable(
ewmaar_monte_carlo_resumo %>%
mutate(across(where(is.numeric), \(x) round(x, digits = 4))),
caption = "Pontos fora de controle por Φ₁ e λ",
colnames = c("Valores de λ", "Valores de Φ₁", "Média", "Mínimo", "Máximo")
)Pela análise gráfica das médias de FPFCs (Fração de Pontos Fora de Controle) notamos que o \(\lambda\) encontrado por otimização ultrapassa o limite de 5% de pontos fora de controle para \(\Phi_1 = 0.2\).
E, de uma forma geral, possui um resultado menos satisfatório que o EWMA sobre os resíduos.
ggplotly(
ewmaar_monte_carlo_resumo %>%
ggplot(aes(
x = h1_phi, y = mean, color = lambda
)) +
geom_line() +
geom_point() +
labs(
x = "Valores de Φ",
y = "Fração de pontos fora de controle",
color = "Valores de λ",
title = "FPFCs por Φ e λ"
) +
theme.base
) %>%
plotly.baseE, a distribuição dos PFCs.
ggplotly(
ewmaar_monte_carlo %>%
ggplot(aes(
x = factor(h1_phi),
y = fracao_fora_de_controle,
fill = lambda
)) +
geom_boxplot() +
labs(
x = "Valores de Φ",
y = "Fração de pontos fora de controle",
fill = "Valores de λ",
title = "FPFCs por Φ e λ"
) +
theme.base
) %>%
plotly.base %>%
layout(boxmode = 'group')O CUSUM (Cumulative Sum) é um método de controle de processos que utiliza a soma acumulada dos resíduos para detectar mudanças no processo.
Neste método são definidas duas constantes, \(H\) e \(K\), que representam o tamanho da mudança que se deseja detectar e a sensibilidade do método, respectivamente. As estatísticas de controle, sendo duas, são definidas como:
\[ \begin{matrix} \text{C}^+_i & = & \max\left[0, x_i - (\mu_0 + K) + C^+_{i-1}\right] \\ \text{C}^-_i & = & \max\left[0, (\mu_0 - K) - x_i + C^-_{i-1}\right] \\ \end{matrix} \]
onde \(x_i\) é a observação no instante \(i\), \(\mu_0\) é a média do processo, \(C^+_0 = C^-_0 = 0\) e \(i = 1, 2, \ldots, n\).
A constante \(K\), chamada de valor de referência, geralmente, segundo Montgomery (2009), é definida como a metade entre o \(\mu_0\) alvo e o valor fora de controle de média \(\mu_1\) que queremos detectar.
Já a constante \(H\), chamada de limite de decisão, é definida como o valor que a estatística de controle deve atingir para podermos detectar a mudança no processo. Para Montgomery (2009), um valor razoável para \(H\) é \(5\sigma\). Assim, se \(\text{C}^+_i \geq H\) ou \(\text{C}^-_i \leq -H\), então o processo está fora de controle.
cusum_qcc <- function(amostra_inicial_,
desvio_detectavel = 1,
intervalo_de_decisao = 5,
amostra_subsequente = NULL) {
arguments <- list(
data = amostra_inicial_,
se.shift = desvio_detectavel,
decision.interval = intervalo_de_decisao,
plot = F
)
if (!is.null(amostra_subsequente)) {
arguments$newdata <- amostra_subsequente
}
cu <- do.call(cusum, arguments)
registros <- if (is.null(amostra_subsequente)) {
seq(1, nH0)
} else {
seq(nH0 + 1, nH0 + nH1)
}
pos <- cu[["pos"]][registros]
neg <- cu[["neg"]][registros]
fora_de_controle_pos <- pos > intervalo_de_decisao
fora_de_controle_neg <- neg < -intervalo_de_decisao
total_fora_de_controle <- sum(ifelse(
fora_de_controle_pos,
fora_de_controle_pos,
fora_de_controle_neg
))
fracao_fora_de_controle <- total_fora_de_controle / length(registros)
list(
pos = pos,
neg = neg,
fora_de_controle_pos = fora_de_controle_pos,
fora_de_controle_neg = fora_de_controle_neg,
total_fora_de_controle = total_fora_de_controle,
fracao_fora_de_controle = fracao_fora_de_controle
)
}KQueremos que, sob \(H_1\), nas 100
amostras iniciais a probabilidade de um ponto fora de controle seja de,
nominalmente, 5%. Portanto, vamos buscar um valor para K
onde isso acontece.
desvio_detectavel_otimo_optimize <- function(amostra_inicial) {
for (x in seq(0, 3, by = 0.1)) {
fora <- cusum_qcc(amostra_inicial, desvio_detectavel = x)$total_fora_de_controle
if (fora == 5) {
return(x)
}
}
return(NA)
}
cusum_k_otimo <- cache_dados("cusum-k-otimo", \() {
data.frame(id = 1:1000) %>%
rowwise() %>%
mutate(
h0_amostra = list(barma.sim(
n = nH0, phi = phi_parametro, seed = id
)),
h0_phi = list(barma.phi_estimado(h0_amostra)),
h0_residuos = list(barma.residuos(h0_amostra, h0_phi)),
minimo = desvio_detectavel_otimo_optimize(h0_residuos)
)
})ggplotly(
cusum_k_otimo %>%
filter(!is.na(minimo)) %>%
ggplot(aes(
x = 0,
y = minimo,
group = 0,
fill = "red"
)) +
geom_boxplot() +
labs(x = element_blank(), y = "Valor ótimo de K", title = "Valor ótimo de K por simulação") +
theme.base +
theme.no_legend
) %>%
plotly.basevalor_otimo_k <- mean(cusum_k_otimo$minimo, na.rm = TRUE)
kable(
cusum_k_otimo %>%
filter(!is.na(minimo)) %>%
group_by() %>%
summarise(
"Média" = mean(minimo),
"Desvio padrão" = sd(minimo),
"Amostras" = n(),
),
caption = "Valor ótimo para K",
align = 'c',
digits = 3
)| Média | Desvio padrão | Amostras |
|---|---|---|
| 0.507 | 0.193 | 260 |
cumsum_k_monte_carlo <- cache_dados("simulacao-cusum-k", function() {
gerador_monte_carlo(parametros = list(desvio_detectavel = seq(0.1, 1.0, by = 0.1))) %>%
mutate(
qcc = list(
cusum_qcc(
h0_residuos,
desvio_detectavel = desvio_detectavel,
amostra_subsequente = h1_residuos
)
),
pos = list(qcc$pos),
neg = list(qcc$neg),
fora_de_controle_pos = list(qcc$fora_de_controle_pos),
fora_de_controle_neg = list(qcc$fora_de_controle_neg),
total_fora_de_controle = qcc$total_fora_de_controle,
fracao_fora_de_controle = qcc$fracao_fora_de_controle
) %>%
select(-qcc)
})K ótimoVamos analisar o CUSUM para o valor ótimo de K
encontrado.
cusum_para_k_otimo <- cumsum_k_monte_carlo %>%
filter(k == 1 & desvio_detectavel == round(valor_otimo_k, 1)) %>%
rename(`Φ₁` = h1_phi) %>%
mutate(observacao = list(1:nH1))
ggplotly(
cusum_para_k_otimo %>%
unnest(
cols = c(
`Φ₁`,
observacao,
pos,
neg,
h1_residuos,
fora_de_controle_pos,
fora_de_controle_neg
)
) %>%
ggplot() +
geom_hline(aes(yintercept = 5, color = "Limite de decisão"), linetype = "dashed") +
geom_hline(aes(
yintercept = -5, color = "Limite de decisão"
), linetype = "dashed") +
geom_line(aes(
x = observacao, y = pos, color = "N+"
)) +
geom_line(aes(
x = observacao, y = neg, color = "N-"
)) +
geom_point(
data = . %>% filter(fora_de_controle_pos),
aes(x = observacao, y = pos, color = "Fora de controle"),
size = 1,
shape = 4
) +
geom_point(
data = . %>% filter(fora_de_controle_neg),
aes(x = observacao, y = neg, color = "Fora de controle"),
size = 1,
shape = 4
) +
labs(x = "Observação", y = "Caracterísitca", color = "Legenda") +
scale_color_manual(
values = c(
"N+" = adjustcolor("blue", alpha.f = 0.6),
"N-" = adjustcolor("darkgreen", alpha.f = 0.6),
"Fora de controle" = adjustcolor("#f42f2f", alpha.f = 0.6),
"Limite de decisão" = adjustcolor("red", alpha.f = 0.5)
)
) +
labs(title = "CUSUM para resíduos βAR(1)", ) +
theme.base +
theme(legend.position = "bottom", legend.title = element_blank()) +
facet_wrap(vars(`Φ₁`), ncol = 2, labeller = "label_both")
) %>%
plotly.baseEm média, diminuir o K aumenta a sensibilidade do CUSUM
para detectar mudanças no processo.
cumsum_k_monte_carlo_resumo <- cumsum_k_monte_carlo %>%
group_by(desvio_detectavel, h1_phi) %>%
summarise(
mean = mean(fracao_fora_de_controle),
min = min(fracao_fora_de_controle),
max = max(fracao_fora_de_controle),
.groups = "drop"
)
datatable(
cumsum_k_monte_carlo_resumo %>%
mutate(across(where(is.numeric), \(x) round(x, digits = 4))),
caption = "Pontos fora de controle por Φ₁ e K, com H = 5",
colnames = c("Valores de K", "Valores de Φ₁", "Média", "Mínimo", "Máximo")
)Análise gráfica das médias de FPFCs.
ggplotly(
cumsum_k_monte_carlo_resumo %>%
ggplot(aes(
x = h1_phi,
y = mean,
color = factor(desvio_detectavel)
)) +
geom_line() +
geom_point() +
labs(
x = "Valores de Φ₁",
y = "Fração de pontos fora de controle",
color = "Valores de K",
title = "FPFCs por Φ₁ e K com H = 5"
) +
theme.base
) %>%
plotly.baseE, a distribuição dos PFCs.
ggplotly(
cumsum_k_monte_carlo %>%
ggplot(aes(
x = factor(h1_phi),
y = fracao_fora_de_controle,
fill = factor(desvio_detectavel)
)) +
geom_boxplot() +
labs(
x = "Valores de Φ₁",
y = "Fração de pontos fora de controle",
fill = "Valores de K",
title = "FPFCs por Φ₁ e K com H = 5"
) +
theme.base
) %>%
plotly.base %>%
layout(boxmode = 'group')HSemelhante ao valor ótimo para K, queremos que, sob
\(H_1\), nas 100 amostras iniciais a
probabilidade de um ponto fora de controle seja de, nominalmente, 5%.
Portanto, vamos buscar um valor para H onde isso
acontece.
intervalo_de_decisao_otimo_optimize <- function(amostra_inicial) {
for (x in seq(1, 6, by = 1)) {
fora <- cusum_qcc(amostra_inicial, intervalo_de_decisao = x)$total_fora_de_controle
if (fora == 5) {
return(x)
}
}
return(NA)
}
cusum_h_otimo <- cache_dados("cusum-h-otimo", \() {
data.frame(id = 1:1000) %>%
rowwise() %>%
mutate(
h0_amostra = list(barma.sim(
n = nH0, phi = phi_parametro, seed = id
)),
h0_phi = list(barma.phi_estimado(h0_amostra)),
h0_residuos = list(barma.residuos(h0_amostra, h0_phi)),
minimo = intervalo_de_decisao_otimo_optimize(h0_residuos)
)
})ggplotly(
cusum_h_otimo %>%
filter(!is.na(minimo)) %>%
ggplot(aes(
x = 0,
y = minimo,
group = 0,
fill = "red"
)) +
geom_boxplot() +
labs(x = element_blank(), y = "Valor ótimo de K", title = "Valor ótimo de K por simulação") +
theme.base +
theme.no_legend
) %>%
plotly.basevalor_otimo_h <- mean(cusum_h_otimo$minimo, na.rm = TRUE)
kable(
cusum_h_otimo %>%
filter(!is.na(minimo)) %>%
group_by() %>%
summarise(
"Média" = mean(minimo),
"Desvio padrão" = sd(minimo),
"Amostras" = n(),
),
caption = "Valor ótimo para H",
align = 'c',
digits = 3
)| Média | Desvio padrão | Amostras |
|---|---|---|
| 3.16 | 0.712 | 125 |
cumsum_h_monte_carlo <- cache_dados("simulacao-cusum-h", function() {
gerador_monte_carlo(parametros = list(intervalo_de_decisao = seq(1, 6, by = 1))) %>%
mutate(
qcc = list(
cusum_qcc(
h0_residuos,
intervalo_de_decisao = intervalo_de_decisao,
amostra_subsequente = h1_residuos
)
),
pos = list(qcc$pos),
neg = list(qcc$neg),
fora_de_controle_pos = list(qcc$fora_de_controle_pos),
fora_de_controle_neg = list(qcc$fora_de_controle_neg),
total_fora_de_controle = qcc$total_fora_de_controle,
fracao_fora_de_controle = qcc$fracao_fora_de_controle
) %>%
select(-qcc)
})K ótimoVamos analisar o CUSUM para o valor ótimo de K
encontrado.
cusum_para_h_otimo <- cumsum_h_monte_carlo %>%
filter(k == 1 & intervalo_de_decisao == round(valor_otimo_h)) %>%
rename(`Φ₁` = h1_phi) %>%
mutate(observacao = list(1:nH1))
ggplotly(
cusum_para_h_otimo %>%
unnest(
cols = c(
`Φ₁`,
observacao,
pos,
neg,
h1_residuos,
fora_de_controle_pos,
fora_de_controle_neg
)
) %>%
ggplot() +
geom_hline(aes(yintercept = 5, color = "Limite de decisão"), linetype = "dashed") +
geom_hline(aes(
yintercept = -5, color = "Limite de decisão"
), linetype = "dashed") +
geom_line(aes(
x = observacao, y = pos, color = "N+"
)) +
geom_line(aes(
x = observacao, y = neg, color = "N-"
)) +
geom_point(
data = . %>% filter(fora_de_controle_pos),
aes(x = observacao, y = pos, color = "Fora de controle"),
size = 1,
shape = 4
) +
geom_point(
data = . %>% filter(fora_de_controle_neg),
aes(x = observacao, y = neg, color = "Fora de controle"),
size = 1,
shape = 4
) +
labs(x = "Observação", y = "Caracterísitca", color = "Legenda") +
scale_color_manual(
values = c(
"N+" = adjustcolor("blue", alpha.f = 0.6),
"N-" = adjustcolor("darkgreen", alpha.f = 0.6),
"Fora de controle" = adjustcolor("#f42f2f", alpha.f = 0.8),
"Limite de decisão" = adjustcolor("red", alpha.f = 0.5)
)
) +
labs(title = "CUSUM para resíduos βAR(1)", ) +
theme.base +
theme(legend.position = "bottom", legend.title = element_blank()) +
facet_wrap(vars(`Φ₁`), ncol = 2, labeller = "label_both")
) %>%
plotly.baseEm média, diminuir o H aumenta a sensibilidade do CUSUM
para detectar mudanças no processo.
cumsum_h_monte_carlo_resumo <- cumsum_h_monte_carlo %>%
group_by(intervalo_de_decisao, h1_phi) %>%
summarise(
mean = mean(fracao_fora_de_controle),
min = min(fracao_fora_de_controle),
max = max(fracao_fora_de_controle),
.groups = "drop"
)
datatable(
cumsum_h_monte_carlo_resumo %>%
mutate(across(where(is.numeric), \(x) round(x, digits = 4))),
caption = "Pontos fora de controle por Φ₁ e H, com K = 1",
colnames = c("Valores de H", "Valores de Φ₁", "Média", "Mínimo", "Máximo")
)Análise gráfica das médias de FPFCs.
ggplotly(
cumsum_h_monte_carlo_resumo %>%
ggplot(aes(
x = h1_phi,
y = mean,
color = factor(intervalo_de_decisao)
)) +
geom_line() +
geom_point() +
labs(
x = "Valores de Φ₁",
y = "Fração de pontos fora de controle",
color = "Valores de H",
title = "FPFCs por Φ₁ e H, com K = 1",
) +
theme.base
) %>%
plotly.baseE, a distribuição dos PFCs.
ggplotly(
cumsum_h_monte_carlo %>%
ggplot(aes(
x = factor(h1_phi),
y = fracao_fora_de_controle,
fill = factor(intervalo_de_decisao)
)) +
geom_boxplot() +
labs(
x = "Valores de Φ₁",
y = "Fração de pontos fora de controle",
fill = "Valores de K",
title = "FPFCs por Φ₁ e K com H = 5"
) +
theme.base
) %>%
plotly.base %>%
layout(boxmode = 'group')Vamos verificar se existe alguma combinação ótima entre os valores de
K e H para detectar pontos fora de
controle.
comb_cusum <- cache_dados("simulacao-cusum-k-h", function() {
gerador_monte_carlo(parametros = expand.grid(
desvio_detectavel = seq(0.6, 1.8, by = 0.2),
intervalo_de_decisao = seq(3, 6, by = 1)
)) %>%
mutate(
qcc = list(
cusum_qcc(
h0_residuos,
desvio_detectavel = desvio_detectavel,
intervalo_de_decisao = intervalo_de_decisao,
amostra_subsequente = h1_residuos
)
),
pos = list(qcc$pos),
neg = list(qcc$neg),
fora_de_controle_pos = list(qcc$fora_de_controle_pos),
fora_de_controle_neg = list(qcc$fora_de_controle_neg),
total_fora_de_controle = qcc$total_fora_de_controle,
fracao_fora_de_controle = qcc$fracao_fora_de_controle
) %>%
select(-qcc)
})comb_cusum_resumo <- comb_cusum %>%
group_by(desvio_detectavel, intervalo_de_decisao, h1_phi) %>%
summarise(
mean = mean(fracao_fora_de_controle),
min = min(fracao_fora_de_controle),
max = max(fracao_fora_de_controle),
.groups = "drop"
) %>%
mutate(
desvio_detectavel = as.factor(desvio_detectavel),
intervalo_de_decisao = as.factor(intervalo_de_decisao),
)
datatable(
comb_cusum_resumo %>%
mutate(across(where(is.numeric), \(x) round(x, digits = 4))),
caption = "Pontos fora de controle",
colnames = c(
"Valores de K",
"Valores de H",
"Valores de Φ₁",
"Média",
"Mínimo",
"Máximo"
)
)ggplotly(
comb_cusum_resumo %>%
group_by(desvio_detectavel, intervalo_de_decisao) %>%
summarise(
h1_phi = list(h1_phi),
mean = list(mean),
.groups = "drop"
) %>%
filter(map_lgl(mean, ~ .x[1] >= 0.05 & .x[1] < 0.06)) %>%
unnest(cols = c(h1_phi, mean)) %>%
ggplot(
aes(
x = h1_phi,
y = mean,
color = desvio_detectavel,
linetype = intervalo_de_decisao
)
) +
geom_line() +
geom_point() +
labs(
x = "Valores de Φ₁",
y = "Fração de pontos fora de controle",
color = "Valores de K",
linetype = "Valores de H",
title = "FPFCs por Φ₁, K e H"
) +
theme.base
) %>%
plotly.baseDe uma forma geral, observa-se que aumentando K e
H diminui a sensibilidade do CUSUM para detectar pontos
fora de controle, e observa-se um leve aumento no poder quando
combinamos H = 5 e K = 0.8, com um leve
aumento na variabilidade.
Vamos comparar os algoritmos EWMA e CUSUM(K) e CUSUM(H) para diferentes valores de \(\Phi_1\).
Aqui chamamos de CUSUM(K) as estatísticas onde o \(H\) for mantido constante e o \(K\) é variado, e o inverso para CUSUM(H).
estatisticas_resumo_combinadas <- bind_rows(
cumsum_k_monte_carlo_resumo %>%
filter(desvio_detectavel == 0.9) %>%
mutate(
desvio_detectavel = as.factor(desvio_detectavel),
algoritmo = "CUSUM(K)"
) %>%
rename(parametro = desvio_detectavel),
cumsum_h_monte_carlo_resumo %>%
filter(intervalo_de_decisao == 4) %>%
mutate(
intervalo_de_decisao = as.factor(intervalo_de_decisao),
algoritmo = "CUSUM(H)"
) %>%
rename(parametro = intervalo_de_decisao),
comb_cusum_resumo %>%
filter(desvio_detectavel == 0.8 & intervalo_de_decisao == 5) %>%
mutate(`(K,H)` = as.factor(
paste(intervalo_de_decisao, desvio_detectavel, sep = ";")
), algoritmo = "CUSUM(K,H)") %>%
rename(parametro = `(K,H)`),
ewma_monte_carlo_resumo %>%
filter(lambda == 0.7) %>%
rename(parametro = lambda) %>% mutate(algoritmo = "EWMA(λ)")
) %>%
select(algoritmo, parametro, h1_phi, mean) %>%
mutate(parametro = as.factor(parametro))ggplotly(
estatisticas_resumo_combinadas %>%
ggplot() +
geom_line(aes(
x = h1_phi,
y = mean,
color = parametro,
linetype = algoritmo
)) +
geom_point(aes(
x = h1_phi,
y = mean,
color = parametro,
shape = algoritmo
)) +
labs(
x = "Valores de Φ₁",
y = "Fração de pontos fora de controle",
linetype = "Algoritmo",
title = "Fração de pontos fora de controle",
color = "Valores dos parâmetros",
shape = "Algoritmo"
) +
theme.base
) %>%
plotly.baseestatisticas_combinadas <- bind_rows(
cumsum_k_monte_carlo %>%
filter(desvio_detectavel == 0.9) %>%
mutate(
desvio_detectavel = as.factor(desvio_detectavel),
algoritmo = "CUSUM(K)"
) %>%
rename(parametro = desvio_detectavel),
cumsum_h_monte_carlo %>%
filter(intervalo_de_decisao == 4) %>%
mutate(
intervalo_de_decisao = as.factor(intervalo_de_decisao),
algoritmo = "CUSUM(H)"
) %>%
rename(parametro = intervalo_de_decisao),
comb_cusum %>%
filter(desvio_detectavel == 0.8 & intervalo_de_decisao == 5) %>%
mutate(`(K,H)` = as.factor(
paste(intervalo_de_decisao, desvio_detectavel, sep = ";")
), algoritmo = "CUSUM(K,H)") %>%
rename(parametro = `(K,H)`),
ewma_monte_carlo %>%
mutate(lambda = as.factor(lambda), algoritmo = "EWMA(λ)") %>%
filter(lambda == 0.7) %>%
rename(parametro = lambda)
) %>%
select(algoritmo, parametro, h1_phi, fracao_fora_de_controle) %>%
mutate(parametro = as.factor(parametro), h1_phi = as.factor(h1_phi))ggplotly(
estatisticas_combinadas %>%
ggplot() +
geom_boxplot(
aes(
x = h1_phi,
y = fracao_fora_de_controle,
fill = parametro,
shape = algoritmo
),
) +
labs(
x = "Valores de Φ₁",
y = "Fração de pontos fora de controle",
fill = "Valores dos parâmetros",
linetype = "Algoritmo",
title = "Fração de pontos fora de controle",
color = "Valores dos parâmetros"
) +
facet_wrap( ~ algoritmo) +
theme.base
) %>%
plotly.base %>%
layout(boxmode = 'group')A partir da análise da média dos pontos fora de controle, observa-se que o CUSUM(K), em média, possui poder semelhante ao CUSUM(H), enquanto o EWMA é inferior aos dois.
Já a análise gráfica dos boxplots mostra que o CUSUM(K) possui uma variabilidade, em geral, maior que o CUSUM(H), enquanto o EWMA possui, em geral, a menor variabilidade entre os algoritmos.
Vamos combinar EWMA e CUSUM e ver o que acontece.
comb_monte_carlo <- cache_dados("simulacao-ewma-cusum", function() {
gerador_monte_carlo(
parametros = expand.grid(
desvio_detectavel = seq(0.2, 1.6, by = 0.2),
intervalo_de_decisao = seq(1, 6, by = 1),
lambda = seq(0.1, 1.0, by = 0.1)
),
numero_de_execucoes = 50
) %>%
mutate(
"(K;H;λ)" = factor(
paste(desvio_detectavel, intervalo_de_decisao, lambda, sep = ";")
),
ewma = list(ewma_qcc(h0_residuos, lambda, h1_residuos)),
cumsum = list(
cusum_qcc(
h0_residuos,
desvio_detectavel = desvio_detectavel,
intervalo_de_decisao = intervalo_de_decisao,
amostra_subsequente = h1_residuos
)
),
fora_de_controle_coop = list(
ewma$fora_de_controle |
(cumsum$fora_de_controle_pos | cumsum$fora_de_controle_neg)
),
fracao_fora_de_controle_coop = sum(fora_de_controle_coop) / nH1,
fora_de_controle_opos = list(
ewma$fora_de_controle &
(cumsum$fora_de_controle_pos | cumsum$fora_de_controle_neg)
),
fracao_fora_de_controle_opos = sum(fora_de_controle_opos) / nH1
) %>%
select(-ewma, -cumsum)
})Vamos verificar como se comporta a combinação dos algoritmos EWMA e CUSUM, ou seja, se um dos algoritmos detectar um ponto fora de controle, então o processo é considerado fora de controle.
Nesta forma mantemos o nível de significância para ambos os algoritmos em 5%.
comb_coop_monte_carlo_resumo <- comb_monte_carlo %>%
group_by(`(K;H;λ)`, h1_phi, desvio_detectavel, lambda) %>%
summarise(
mean = mean(fracao_fora_de_controle_coop),
min = min(fracao_fora_de_controle_coop),
max = max(fracao_fora_de_controle_coop),
.groups = "drop"
) %>%
group_by(`(K;H;λ)`) %>%
summarise(
h1_phi = list(h1_phi),
desvio_detectavel = list(desvio_detectavel),
lambda = list(lambda),
mean = list(mean),
min = list(min),
max = list(max),
.groups = "drop"
) %>%
filter(map_lgl(mean, ~ .x[1] >= 0.05 & .x[1] < 0.06)) %>%
unnest(cols = c(h1_phi, desvio_detectavel, lambda, mean, min, max))
datatable(
comb_coop_monte_carlo_resumo %>%
select(-`(K;H;λ)`) %>%
mutate(across(where(is.numeric), \(x) round(x, digits = 4))),
caption = "Pontos fora de controle por Φ₁ e (K;H;λ)",
colnames = c("Valores de Φ₁", "K", "H", "λ", "Média", "Mínimo", "Máximo")
)ggplotly(
comb_coop_monte_carlo_resumo %>%
group_by(`(K;H;λ)`) %>%
summarise(
h1_phi = list(h1_phi),
mean = list(mean),
.groups = "drop"
) %>%
filter(map_lgl(mean, ~ .x[1] >= 0.05 & .x[1] < 0.06)) %>%
unnest(cols = c(h1_phi, mean)) %>%
ggplot(aes(
x = h1_phi, y = mean, color = `(K;H;λ)`
)) +
geom_line() +
geom_point() +
labs(
x = "Valores de Φ₁",
y = "Fração de pontos fora de controle",
color = "Constantes (K;H;λ)",
title = "FPFC por Φ₁ e (K;H;λ) (coop)"
) +
theme.base
) %>%
plotly.baseVamos verificar como se comporta a combinação dos algoritmos EWMA vs. CUSUM, ou seja, apenas se ambos os algoritmos detectarem um ponto fora de controle, então o processo é considerado fora de controle.
Fazendo isso, diminui-se o nível de significância para cada um dos algoritmos. Ainda queremos que o nível de significância final seja de 5%, para tanto buscamos uma relação \(0.05 = 1 - (1 - \alpha_{EWMA})(1 - \alpha_{CUSUM})\), obtida a partir da significância para testes simultâneos.
comb_opos_monte_carlo_resumo <- comb_monte_carlo %>%
group_by(`(K;H;λ)`, h1_phi, desvio_detectavel, lambda) %>%
summarise(
mean = mean(fracao_fora_de_controle_opos),
min = min(fracao_fora_de_controle_opos),
max = max(fracao_fora_de_controle_opos),
.groups = "drop"
) %>%
group_by(`(K;H;λ)`) %>%
summarise(
h1_phi = list(h1_phi),
desvio_detectavel = list(desvio_detectavel),
lambda = list(lambda),
mean = list(mean),
min = list(min),
max = list(max),
.groups = "drop"
) %>%
filter(map_lgl(mean, ~ .x[1] >= 0.05 & .x[1] < 0.06)) %>%
unnest(cols = c(h1_phi, desvio_detectavel, lambda, mean, min, max))
datatable(
comb_opos_monte_carlo_resumo %>%
select(-`(K;H;λ)`) %>%
mutate(across(where(is.numeric), \(x) round(x, digits = 4))),
caption = "Pontos fora de controle por Φ₁ e (K;H;λ)",
colnames = c("Valores de Φ₁", "K", "H", "λ", "Média", "Mínimo", "Máximo")
)ggplotly(
comb_opos_monte_carlo_resumo %>%
group_by(`(K;H;λ)`) %>%
summarise(
h1_phi = list(h1_phi),
mean = list(mean),
.groups = "drop"
) %>%
filter(map_lgl(mean, ~ .x[1] >= 0.05 & .x[1] < 0.052)) %>%
unnest(cols = c(h1_phi, mean)) %>%
ggplot(aes(
x = h1_phi, y = mean, color = `(K;H;λ)`
)) +
geom_line() +
geom_point() +
labs(
x = "Valores de Φ₁",
y = "Fração de pontos fora de controle",
color = "Constantes (K;H;λ)",
title = "FPFC por Φ₁ e (K;H;λ) (opos)"
) +
theme.base
) %>%
plotly.baseestatisticas_combinadas_resumo_ew_sum <- bind_rows(
estatisticas_resumo_combinadas,
comb_coop_monte_carlo_resumo %>%
filter(`(K;H;λ)` == "1.4;4;0.8") %>%
mutate(parametro = as.factor(`(K;H;λ)`)) %>%
mutate(algoritmo = "EW-SUM(coop)") %>%
select(algoritmo, parametro, h1_phi, mean),
comb_opos_monte_carlo_resumo %>%
filter(`(K;H;λ)` == "1;4;0.1") %>%
mutate(parametro = as.factor(`(K;H;λ)`)) %>%
mutate(algoritmo = "EW-SUM(opos)") %>%
select(algoritmo, parametro, h1_phi, mean)
)ggplotly(
estatisticas_combinadas_resumo_ew_sum %>%
ggplot() +
geom_line(aes(
x = h1_phi,
y = mean,
color = parametro,
linetype = algoritmo
)) +
geom_point(aes(
x = h1_phi,
y = mean,
color = parametro,
shape = algoritmo
)) +
labs(
x = "Valores de Φ₁",
y = "Fração de pontos fora de controle",
linetype = "Algoritmo",
title = "Fração de pontos fora de controle",
color = "Valores dos parâmetros",
shape = "Algoritmo"
) +
theme.base
) %>%
plotly.base